function crit=estimMarriage20191030(thetaR)
tic;
global beta  lambdaH lambdaHNI A delta sigma AgeDiffCoeff OutsideCoeff Glob THETA Param Nationality gamma0 gamma1
global lambdaNatMen lambdaEducF lambdaAgeM lambdaNation lambdaNorth lambdaLIM2 lambdaGIM lambdaGIM2 lambdaLIW lambdaLIW2 lambdaGIW lambdaGIW2
global indactive indA indbeta indlambda indlambdaH indsigma indlambdaH2 indOutside inddelta indgamma
global alpha AgeParam  NationParam indAgeParam indNation stderr
global dimS dimProv dimNum dimAge dimEduc dimNation dimHouse ageGrid educGrid typicalSurplus NorthSouth
global Match1 Match2 SinglesM1 SinglesF1 SinglesM2 SinglesF2 Weight ObsFertility
global PredictMarriage1 PredictMarriage2 PredictSinglesM1 PredictSinglesM2 PredictSinglesF1 PredictSinglesF2 ProvDistribMen ProvDistribWomen
global mincrit first 
global Italian ItalianLowLow ItalianLowHigh ItalianHighLow ItalianHighHigh  ItalianMen ItalianWomen NonItalian Nationality IndMarriage
global PredictMarriage SURPLUS CoeffDivorce gamma VecTotSingles New SimFertility
theta=THETA;
theta(indactive)=thetaR;
beta=theta(indbeta);    
sigma=abs(theta(indsigma));
lambdaNatMen=abs(theta(indlambda(1)));
lambdaAgeM=theta(indlambda(2));
lambdaEducF=theta(indlambda(3));
lambdaNation=theta(indlambda(4:10));
lambdaNorth=theta(indlambda(11));
delta=abs(theta(inddelta));
lambdaH=theta(indlambdaH(1:3));
lambdaHNI=abs(theta(indlambdaH2(1:3)));

A=theta(indA(1:4));
A=reshape(A,2,2);
OutsideCoeff=zeros(dimProv,2);
OutsideCoeff(1:dimProv,1)=[theta(indOutside(1));theta(indOutside(1))*theta(indOutside(2:dimProv))];
OutsideCoeff(1:dimProv,2)=[theta(indOutside(1));theta(indOutside(1))*theta(indOutside(2:dimProv))]*theta(indOutside(dimProv+1));
AgeParam=[1;theta(indAgeParam)];
AgeParam=reshape(AgeParam,dimAge,dimAge);

alphaOffDiag=constrain(theta(indNation(1:28)),0,1).^lambdaNation(3);

Nationality=zeros(dimNation,dimNation);
Nationality(logical(triu(ones(dimNation,dimNation),1)))=alphaOffDiag;
Nationality=Nationality'+Nationality;
Nationality(logical(eye(dimNation)))=1;
Nationality(1,3:end)=Nationality(1,3:end).^lambdaNatMen;
gamma0=theta(indgamma(1));
gamma1=theta(indgamma(2));

Param.A=A;
Param.AgeParam=AgeParam;
Param.beta=beta;
Param.lambdaNatMen=lambdaNatMen;
Param.lambdaAgeM=lambdaAgeM;
Param.lambdaEducF=lambdaEducF;
Param.delta=delta;
Param.sigma=sigma;
Param.AgeDiff=AgeDiffCoeff;
Param.Outside=OutsideCoeff;
Param.OutsideNI=[0;theta(indNation(29:end))];
Param.alpha=alpha;
Param.lambdaH=lambdaH;
Param.lambdaHNI=lambdaHNI;
Param.Nationality=Nationality;
Param.lambdaNation=lambdaNation;
Param.lambdaNorth=lambdaNorth;

dimProv=Glob.dimProv;
dimAge=Glob.dimAge;
dimEduc=Glob.dimEduc;
dimNation=Glob.dimNation;
DefaultCoeff=Glob.DefaultCoeff;


PredictMarriage=zeros(dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimProv*dimNum,2);
STORE=zeros(dimS,14,dimProv*dimNum);
STORE_=zeros(dimS,14,dimProv*dimNum);
STOREAfter=zeros(dimS,14,dimProv*dimNum);
STOREAfter_=zeros(dimS,14,dimProv*dimNum);
STORE3YB=zeros(dimS,8,dimProv*dimNum);
STORE3YB_=zeros(dimS,8,dimProv*dimNum);
STORE3YA=zeros(dimS,8,dimProv*dimNum);
STORE3YA_=zeros(dimS,8,dimProv*dimNum);


PROV=kron(1:dimProv,ones(1,dimNum));
NUM=kron(ones(1,dimProv),1:dimNum);
gamma=zeros(dimProv*dimNum,1);
Period=1;    

    parfor loop=1:(dimProv*dimNum)
    iprov=PROV(loop);
    num=NUM(loop);
    [Pred_Marriage,TOTAL_SURPLUS,allstore,gamma_]=solveMatch20191030(Param,Glob,iprov,Period,num);
    PredictMarriage(:,:,loop,Period)=Pred_Marriage;
     STORE(:,:,loop)=allstore.STORE;
     STORE_(:,:,loop)=allstore.STORE_;
     STORE3YB(:,:,loop)=allstore.STORE3YB;
     STORE3YB_(:,:,loop)=allstore.STORE3YB_;
     STORE3YA(:,:,loop)=allstore.STORE3YA;
     STORE3YA_(:,:,loop)=allstore.STORE3YA_;   
     gamma(loop)=gamma_;
    end

Period=2;    

    parfor loop=1:(dimProv*dimNum)
    iprov=PROV(loop);
    num=NUM(loop);
    [Pred_Marriage,TOTAL_SURPLUS,allstore,gamma_]=solveMatch20191030(Param,Glob,iprov,Period,num);
    PredictMarriage(:,:,loop,Period)=Pred_Marriage;
    STORE3YC(:,:,loop)=allstore.STORE3YC;
    STORE3YC_(:,:,loop)=allstore.STORE3YC_;   
    STOREAfter(:,:,loop)=allstore.STORE;
    STOREAfter_(:,:,loop)=allstore.STORE_;

    end


PredictMarriage1=PredictMarriage(:,:,:,1);
PredictMarriage1=reshape(PredictMarriage1,dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimNum,dimProv);
PredictMarriage1=permute(PredictMarriage1,[1 2 4 3]);
PredictMarriage1=mean(PredictMarriage1,4); 
PredictMarriage2=PredictMarriage(:,:,:,2);
PredictMarriage2=reshape(PredictMarriage2,dimAge*dimEduc*dimNation,dimAge*dimEduc*dimNation,dimNum,dimProv);
PredictMarriage2=permute(PredictMarriage2,[1 2 4 3]);
PredictMarriage2=mean(PredictMarriage2,4); 

nationGrid=kron(ones(dimAge*dimEduc,1),(1:dimNation)');
ageGrid2=kron(kron(ones(dimEduc,1),ageGrid),ones(dimNation,1));
Italian=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])==1)&(repmat(nationGrid',[dimAge*dimEduc*dimNation 1])==1);
Italian=repmat(Italian,[1 1 dimProv]);
ItalianMenOld=(repmat(nationGrid,[1 dimAge*dimEduc*dimNation])==1)&...
    (repmat(ageGrid2>=37,[1 dimAge*dimEduc*dimNation]))&...
    (repmat(nationGrid',[dimAge*dimEduc*dimNation 1])>2);
ItalianMenOld=repmat(ItalianMenOld,[1 1 dimProv]);

DiffItalian=100*(PredictMarriage1(Italian)-Match1(Italian));DiffItalian=DiffItalian'*DiffItalian;


auxdiff1North=zeros(dimNation,dimNation);
auxdiff1South=zeros(dimNation,dimNation);
for i=1:dimNation
    for j=1:dimNation
        auxInd=IndMarriage(:,:,i,j);
        auxInd=repmat(auxInd,[1 1 dimProv]);
        auxIndNorth=auxInd;
        auxIndSouth=auxInd;
        auxIndSouth(:,:,[1 2 3 4 5 6 11 16 17 18 19 20])=0;
        auxIndNorth(:,:,[7 8 9 10 12 13 14 15 21])=0;
        auxN=100*(sum(PredictMarriage1(auxIndNorth))-sum(Match1(auxIndNorth)));auxdiff1North(i,j)=auxN^2;
        auxS=100*(sum(PredictMarriage1(auxIndSouth))-sum(Match1(auxIndSouth)));auxdiff1South(i,j)=auxS^2;
    end
end
auxdiff1North(1,1)=0;auxdiff1North=sum(auxdiff1North(:));
auxdiff1South(1,1)=0;auxdiff1South=sum(auxdiff1South(:));
DiffItOld=100*(sum(PredictMarriage1(ItalianMenOld))-sum(Match1(ItalianMenOld)));DiffItOld=DiffItOld^2;

DiffMarriage1=DiffItalian+(auxdiff1North+auxdiff1South+DiffItOld);


% Period 2
DiffItalian=100*(PredictMarriage2(Italian)-Match2(Italian));DiffItalian=DiffItalian'*DiffItalian;

auxdiff2North=zeros(dimNation,dimNation);
auxdiff2South=zeros(dimNation,dimNation);
for i=1:dimNation
    for j=1:dimNation
        auxInd=IndMarriage(:,:,i,j);
        auxInd=repmat(auxInd,[1 1 dimProv]);
        auxIndNorth=auxInd;
        auxIndSouth=auxInd;
        auxIndSouth(:,:,[1 2 3 4 5 6 11 16 17 18 19 20])=0;
        auxIndNorth(:,:,[7 8 9 10 12 13 14 15 21])=0;
        auxN=100*(sum(PredictMarriage2(auxIndNorth))-sum(Match2(auxIndNorth)));auxdiff2North(i,j)=auxN^2;
        auxS=100*(sum(PredictMarriage2(auxIndSouth))-sum(Match2(auxIndSouth)));auxdiff2South(i,j)=auxS^2;
    end
end
auxdiff2North(1,1)=0;auxdiff2North=sum(auxdiff2North(:));
auxdiff2South(1,1)=0;auxdiff2South=sum(auxdiff2South(:));



DiffItOld=100*(sum(PredictMarriage2(ItalianMenOld))-sum(Match2(ItalianMenOld)));DiffItOld=DiffItOld^2;
DiffMarriage2=DiffItalian+(auxdiff2North+auxdiff2South+DiffItOld);

% Divorces before the enlargement
auxSTORE3YB=reshape(permute(STORE3YB,[1 3 2]),dimS*dimProv*dimNum,8);
WeightRegB=kron(gamma,ones(dimS,1));
DivorceB=auxSTORE3YB(auxSTORE3YB(:,8)==1,7); % Select those who married 3 years prior
WeightRegB=WeightRegB(auxSTORE3YB(:,8)==1);
auxSTORE3YB_=reshape(permute(STORE3YB_,[1 3 2]),dimS*dimProv*dimNum,8);
WeightRegB_=kron(1-gamma,ones(dimS,1));
DivorceB_=auxSTORE3YB_(auxSTORE3YB_(:,8)==1,7); % Select those who married 3 years prior
WeightRegB_=WeightRegB_(auxSTORE3YB_(:,8)==1);
AllDB=[DivorceB;DivorceB_];%AllD=Divorce;
WeightReg=[WeightRegB;WeightRegB_];
X=[auxSTORE3YB(auxSTORE3YB(:,8)==1,1:6);auxSTORE3YB_(auxSTORE3YB_(:,8)==1,1:6)];

dnat_ita_neweu=((X(:,1)==1)&(X(:,2)==3))|((X(:,1)==3)&(X(:,2)==1));
dnat_ita_others=((X(:,1)==1)&(X(:,2)~=1)&(X(:,2)~=3))|((X(:,2)==1)&(X(:,1)~=1)&(X(:,1)~=3));
dnat_othersnew_othersnew=(X(:,1)>=3)&(X(:,2)>=3);
woman10younger=(X(:,4)<(X(:,3)-10));
woman10older=(X(:,4)>(X(:,3)+10));
sameduc=(X(:,5)==X(:,6));

XX=[dnat_ita_neweu dnat_ita_others dnat_othersnew_othersnew woman10younger woman10older sameduc ones(length(X),1)];


XX(end,std(XX(:,1:6))==0)=1; % To prevent the estimation to crash
CoeffDivorceB = olsmomentW(AllDB,XX,WeightReg);
CoeffDivorceB(isnan(CoeffDivorceB))=DefaultCoeff; 

% Divorces after the enlargement
auxSTORE3YA=reshape(permute(STORE3YA,[1 3 2]),dimS*dimProv*dimNum,8);
WeightRegA=kron(gamma,ones(dimS,1));
DivorceA=auxSTORE3YA(auxSTORE3YA(:,8)==1,7); % Select those who married 3 years prior
WeightRegA=WeightRegA(auxSTORE3YA(:,8)==1);
auxSTORE3YA_=reshape(permute(STORE3YA_,[1 3 2]),dimS*dimProv*dimNum,8);
WeightRegA_=kron(1-gamma,ones(dimS,1));
DivorceA_=auxSTORE3YA_(auxSTORE3YA_(:,8)==1,7); % Select those who married 3 years prior
WeightRegA_=WeightRegA_(auxSTORE3YA_(:,8)==1);
AllDA=[DivorceA;DivorceA_];%AllD=Divorce;
WeightReg=[WeightRegA;WeightRegA_];

X=[auxSTORE3YA(auxSTORE3YA(:,8)==1,1:6);auxSTORE3YA_(auxSTORE3YA_(:,8)==1,1:6)];

dnat_ita_neweu=((X(:,1)==1)&(X(:,2)==3))|((X(:,1)==3)&(X(:,2)==1));
dnat_ita_others=((X(:,1)==1)&(X(:,2)~=1)&(X(:,2)~=3))|((X(:,2)==1)&(X(:,1)~=1)&(X(:,1)~=3));
dnat_othersnew_othersnew=(X(:,1)>=3)&(X(:,2)>=3);
woman10younger=(X(:,4)<(X(:,3)-10));
woman10older=(X(:,4)>(X(:,3)+10));
sameduc=(X(:,5)==X(:,6));

XX=[dnat_ita_neweu dnat_ita_others dnat_othersnew_othersnew woman10younger woman10older sameduc ones(length(X),1)];
XX(1,std(XX(:,1:6))==0)=1; % To prevent the estimation to crash
% CoeffDivorceA = ((XX' * XX) \ XX') * AllDA;
CoeffDivorceA = olsmomentW(AllDA,XX,WeightReg);
CoeffDivorceA(isnan(CoeffDivorceA))=DefaultCoeff; 

% Divorces of couples formed after the enlargement
auxSTORE3YC=reshape(permute(STORE3YC,[1 3 2]),dimS*dimProv*dimNum,8);
WeightRegC=kron(gamma,ones(dimS,1));
DivorceC=auxSTORE3YC(auxSTORE3YC(:,8)==1,7); % Select those who married 3 years prior
WeightRegC=WeightRegC(auxSTORE3YC(:,8)==1);
auxSTORE3YC_=reshape(permute(STORE3YC_,[1 3 2]),dimS*dimProv*dimNum,8);
WeightRegC_=kron(1-gamma,ones(dimS,1));
DivorceC_=auxSTORE3YC_(auxSTORE3YC_(:,8)==1,7); % Select those who married 3 years prior
WeightRegC_=WeightRegC_(auxSTORE3YC_(:,8)==1);
AllDC=[DivorceC;DivorceC_];%AllD=Divorce;
WeightReg=[WeightRegC;WeightRegC_];

X=[auxSTORE3YC(auxSTORE3YC(:,8)==1,1:6);auxSTORE3YC_(auxSTORE3YC_(:,8)==1,1:6)];

dnat_ita_neweu=((X(:,1)==1)&(X(:,2)==3))|((X(:,1)==3)&(X(:,2)==1));
dnat_ita_others=((X(:,1)==1)&(X(:,2)~=1)&(X(:,2)~=3))|((X(:,2)==1)&(X(:,1)~=1)&(X(:,1)~=3));
dnat_othersnew_othersnew=(X(:,1)>=3)&(X(:,2)>=3);
woman10younger=(X(:,4)<(X(:,3)-10));
woman10older=(X(:,4)>(X(:,3)+10));
sameduc=(X(:,5)==X(:,6));

XX=[dnat_ita_neweu dnat_ita_others dnat_othersnew_othersnew woman10younger woman10older sameduc ones(length(X),1)];
XX(1,std(XX(:,1:6))==0)=1; % To prevent the estimation to crash
CoeffDivorceC = olsmomentW(AllDC,XX,WeightReg);
CoeffDivorceC(isnan(CoeffDivorceC))=DefaultCoeff; 


DivorceMoments=Glob.DivorceMoments;
CoeffDivorce=cat(3,CoeffDivorceB,CoeffDivorceA,CoeffDivorceC);

DiffDivorce=(CoeffDivorce-DivorceMoments(:,1,:))./DivorceMoments(:,2,:);
DiffDivorce=DiffDivorce(:)'*DiffDivorce(:)+(CoeffDivorce(1,1,1)>CoeffDivorce(1,1,2))*10+(CoeffDivorce(1,1,2)<0)*20;

%% Fertility
optsurplus1=[squeeze(STORE(:,14,:));squeeze(STORE_(:,14,:))];
optsurplus2=[squeeze(STOREAfter(:,14,:));squeeze(STOREAfter_(:,14,:))];
single1=[squeeze(STORE(:,7,:));squeeze(STORE_(:,7,:))];
single2=[squeeze(STOREAfter(:,7,:));squeeze(STOREAfter_(:,7,:))];

optsurplus1=optsurplus1(not(single1));
optsurplus2=optsurplus2(not(single2));
seed=Glob.seed;
rng(seed,'twister');
ufertility1=randn(length(optsurplus1),1);
ufertility2=randn(length(optsurplus2),1);
fertility1=(gamma0+gamma1*optsurplus1+ufertility1)>0;
fertility2=(gamma0+gamma1*optsurplus2+ufertility2)>0;
Province=repmat(PROV(1:(dimProv*dimNum)),[2*dimS 1]);
NatM1=[squeeze(STORE(:,1,:));squeeze(STORE_(:,1,:))];
NatF1=[squeeze(STORE(:,2,:));squeeze(STORE_(:,2,:))];
NatM2=[squeeze(STOREAfter(:,1,:));squeeze(STOREAfter_(:,1,:))];
NatF2=[squeeze(STOREAfter(:,2,:));squeeze(STOREAfter_(:,2,:))];
NatM1=NatM1(not(single1));
NatF1=NatF1(not(single1));
NatM2=NatM2(not(single2));
NatF2=NatF2(not(single2));
Province1=Province(not(single1));
Province2=Province(not(single2));

EU15=(NatM1==1&NatF1==2);
EU2004=(NatM1==1&NatF1==3);
EUOther=(NatM1==1&NatF1==4);
Africa=(NatM1==1&NatF1==5);
Asia=(NatM1==1&NatF1==6);
America=(NatM1==1&NatF1==7);
OECD=(NatM1==1&NatF1==8);
P1=Province1==1;
for i=2:dimProv-1
 P1=[P1 Province1==i];
end 

X=[EU2004 P1 ones(length(OECD),1)];
Y=fertility1;
anyItalian=logical((NatM1==1|NatF1==1)&(not(NatM1==1&NatF1==1)));
Y=Y(anyItalian);
X=X(anyItalian,:);
b1_ItalMales=inv(X'*X)*X'*Y;

EU15=(NatM1==2&NatF1==1);
EU2004=(NatM1==3&NatF1==1);
EUOther=(NatM1==4&NatF1==1);
Africa=(NatM1==5&NatF1==1);
Asia=(NatM1==6&NatF1==1);
America=(NatM1==7&NatF1==1);
OECD=(NatM1==8&NatF1==1);

X=[EU2004  P1 ones(length(OECD),1)];
Y=fertility1;
anyItalian=logical((NatM1==1|NatF1==1)&(not(NatM1==1&NatF1==1)));
Y=Y(anyItalian);
X=X(anyItalian,:);
b1_ItalFemales=inv(X'*X)*X'*Y;

EU15=(NatM2==1&NatF2==2);
EU2004=(NatM2==1&NatF2==3);
EUOther=(NatM2==1&NatF2==4);
Africa=(NatM2==1&NatF2==5);
Asia=(NatM2==1&NatF2==6);
America=(NatM2==1&NatF2==7);
OECD=(NatM2==1&NatF2==8);
P2=Province2==1;
for i=2:dimProv-1
 P2=[P2 Province2==i];
end 

X=[EU2004 P2 ones(length(OECD),1)];
Y=fertility2;
anyItalian=logical((NatM2==1|NatF2==1)&(not(NatM2==1&NatF2==1)));
Y=Y(anyItalian);
X=X(anyItalian,:);
b2_ItalMales=inv(X'*X)*X'*Y;

EU15=(NatM2==2&NatF2==1);
EU2004=(NatM2==3&NatF2==1);
EUOther=(NatM2==4&NatF2==1);
Africa=(NatM2==5&NatF2==1);
Asia=(NatM2==6&NatF2==1);
America=(NatM2==7&NatF2==1);
OECD=(NatM2==8&NatF2==1);

X=[EU2004 P2 ones(length(OECD),1)];
Y=fertility2;
anyItalian=logical((NatM2==1|NatF2==1)&(not(NatM2==1&NatF2==1)));
Y=Y(anyItalian);
X=X(anyItalian,:);
b2_ItalFemales=inv(X'*X)*X'*Y;
SimFertility=[b1_ItalMales(end) b1_ItalFemales(end) b2_ItalMales(1)-b1_ItalMales(1) b2_ItalFemales(1)-b1_ItalFemales(1)];
DiffFertility=(SimFertility-ObsFertility(1,:))./ObsFertility(2,:);
DiffFertility2=DiffFertility*DiffFertility';

crit=DiffMarriage1+DiffMarriage2+DiffDivorce+DiffFertility2; 

disp(toc)

if stderr==1
    SimM=[PredictMarriage1(:);PredictMarriage2(:);CoeffDivorce(:);SimFertility'];
    auxD=DivorceMoments(:,1,:);
    ObsM=[Match1(:);Match2(:);auxD(:);ObsFertility(1,:)'];
    auxD=DivorceMoments(:,2,:);
    denom=[ones(2*dimAge*dimEduc*dimNation*dimAge*dimEduc*dimNation*dimProv,1)/10;auxD(:);ObsFertility(2,:)'];
    crit=(ObsM-SimM)./denom;
end
    